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We study qubit-mediated energy transfer between two electron reservoirs by adopting a 
numerically-exact influence functional path-integral method. This non-perturbative technique al- 
lows us to study the system's dynamics beyond the weak coupling limit. Our simulations for the 
energy current indicate that perturbative-Markovian Master equation predictions significantly de- 
viate from exact numerical results already at intermediate coupling, npctjji > 0.4, where p is the 
metal (Fermi sea) density of states, taken as a constant, and aijji is the scattering potential energy 
of electrons, between the j and j' states. Markovian Master equation techniques should be there- 
fore used with caution beyond the strictly weak subsystem-bath coupling limit, especially when a 
quantitative knowledge of transport characteristics is desired. 

PACS numbers: 05.60-k, 44.10.+i, 44.40.+a, 73.23.-b 
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I. INTRODUCTION 

Quantum impurity models, including a subsystem in- 
teracting with a reservoir, were proven useful in describ- 
ing and predicting many physical phenomena. The spin- 
boson model representing the dynamics of a single 
charge on two states coupled to a dissipative bath, e.g., 
a solvent, exhibits rich phenomenology, including var- 
ious phase transitions [2;]. It is relevant for modeling 
charge transfer reactions in biological systems Q, pho- 
tosynthesis the Kondo problem for magnetic impuri- 
ties and quantum information processing in super- 
conducting Josephson tunneling junction qubits @ or 
nitrogen- vacancy centers in diamonds @. A variant of 
the spin-boson model is the spin-fermion model, where a 
qubit, referred to as a spin or a two-level system, inter- 
acts with a metallic-fermionic environment. This model 
is also related to the Kondo model [H , only lacking direct 
coupling of the reservoir degrees of freedom to spin-flip 
processes. The generalization of the equilibrium spin- 
fermion model, to include more than one Fermi bath, 
provides a minimal setting for the study of dissipation 
and decoherence effects under the influence of an out-of- 
equilibrium environment ■ 

In this work, we use the two-bath spin-fermion model 
and investigate energy exchange between two metals, me- 
diated by the excitation/relaxation of a nonlinear quan- 
tum system, a qubit. For a scheme of this setup, see 
Fig. Q] Physically, our model can describe the process 
of radiative heat transfer between metals |llT[l3| . and 
it can be realized within a superconducting Josephson 
junction circuit [l(| EH, EH ■ We simulate the energy cur- 
rent characteristics of the nonequilibrium spin-fermion 
model in a large parameter range of coupling strengths by 
means of an influence-functional path-integral (INFPI) 
technique developed in Rcfs. 0, [l3]- This numerically- 
exact method is built about the basic observation, that 
in out-of-equilibrium (and/or finite temperature) situ- 
ations bath correlations have a finite range, allowing 



for their truncation beyond a memory time dictated by 
the voltage-bias and the temperature [1, H, EH . Taking 
advantage of this fact, an iterative-deterministic time- 
evolution scheme can be developed, where convergence 
with respect to the memory length can in principle be 
reached. 

Our main objective here is to explore the qubit- 
mcdiatcd energy current characteristics beyond standard 
perturbative methods. Particularly, we would like to find 
when the Golden-Rule-type Markovian Master equation 
method provides a correct (quantitative or qualitative) 
description of the exact behavior. This task is impor- 
tant since Master equation tools have been extensively 
adopted for studying problems in charge, spin, and en- 
ergy transfer phenomenology in quantum dots and molec- 
ular junctions, see for example [l3l [l9l - l28| . 

The plan of the paper is as follows. In Sec. II we 
present the nonequilibrium spin-fermion model. We pro- 
vide expressions for observables of interest in Sec. III. 
The two methods confronted, INFPI and Markovian 
Master equation, are discussed in Sec. IV, with results 
included in Sec. V. Sec VI. concludes. For simplicity, we 
use the conventions % = 1, electron charge e = 1, and 
Boltzmann constant ks = 1. 



II. MODEL 



The system of interest comprises of two metallic leads, 
v = L, R, prepared at different temperatures but at the 
same chemical potential. These metals are connected 
indirectly, by a nonlinear quantum unit, a two-level sub- 
system. The Hamiltonian includes three contributions, 



H = H S + H F + V, 



(1) 
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FIG. 1. A schematic representation of our model system. 
Electron transfer between the metals is blocked, but energy 
current is flowing through an excitation/de-excitation of the 
intermediate anharmonic (two-state) quantum system. The 
curved arrows represent energy transfer processes between the 
leads and the intermediating subsystem. In our work here we 
set fj,L = PR and take Tl > Tr. 



where 



H F = H L + Hr, H v = e ->' c " -• c '. 



V = Vl + Vr, V„ = cr x ^2 a »,r,v,j' c l,j c ^,j' • ( 2 ) 

3,3' 

The subsystem Hs includes two states, |0) and |1), with 
a tunneling splitting 2A. It can be realized within a non- 
linear resonator mode, or it can represent an impurity 
in a solid-state environment. This subsystem interacts 
with two fermionic reservoirs (Hp) where c£j (c»,j) cre- 
ates (annihilates) an electron at the v = L,R metal lead 
with momentum j, disregarding the electron spin degree 
of freedom. The qubit-metal interaction term V couples 
scattering events within each metal to transitions within 
the subsystem. For simplicity, we assume that the cou- 
pling constants a v are energy independent and real num- 
bers. Note that we do not allow for charge transfer pro- 
cesses between the two metals, assuming the tunneling 
barrier is high. However, energy is transferred between 
the two metals, mediated by the excitation of the inter- 
mediate nonlinear quantum system, see Fig. [TJ Using the 
Hamiltonian form ([2]), the subsystem dynamics and the 
energy current can be readily attained within a Marko- 
vian Master equation, as we explain in Sec. IV. B 

The Hamiltonian ([2]) can be transformed into the stan- 
dard spin-fermion model of zero energy spacing with a 
unitary transformation, 

tfa z U = a x , U*a x U = a z , (3) 

where U — 771 (°x + &z)- The transformed Hamiltonian 
Hsf — U^HU includes a a z -type electron-spin coupling, 



H SF = Aa x + Y £j c l,j c v,: 



v,3 



v,3,3' 



(4) 



In this representation, the dynamics can be conveniently 
simulated with INFPI, a brief discussion is included in 
Sec. IV.A. 



III. OBSERVABLES 

We assume a factorized initial state with the total 
density matrix p(0) = ps(0) ® pf- Here ps denotes 
the reduced density matrix of the subsystem. The 
reservoirs' density matrix at time t = is given by 
Pf = Pl <8> Pr, and these states are canonical, p v 



-P v (H„-y. v N„) 



/Tr„[e 



']. Here we trace over 



the v reservoirs degrees of freedom. In our simulations 
below we take p L — Pr, but assume different initial 
temperatures, T L ^ Tr. We refer to this setup as a 
"noncquilibrium environment" since the two reservoirs 
are prepared in different states. At time t = we put 
into contact the two Fermi baths through the quantum 
subsystem, and follow the evolution of the reduced den- 
sity matrix and the energy current, to steady-state. Since 
the energy current in the system is driven by a tempera- 
ture bias, we can also refer to the energy current here as 
a "heat current". 

The time evolution of the reduced density matrix is ob- 
tained by tracing p over the fermionic reservoirs' degrees 
of freedom, 



p s (t) = Tv F [e- im p(0)e im ] 



(5) 



The definition of the energy current operator is more 
subtle [1^, [3^ , as different-plausible choices provide dis- 
tinct results in the short time limit. At long time, in the 
steady-state limit, these definitions yield the same value. 
Here, we follow the analysis of Ref. [29(, and define the 
energy current operator, e.g., at the left contact, as 



Jv = -[H L -H S ,V L \. 



(6) 



In steady-state the expectation value of the interaction 
is zero, and we reach the relation, 



Tr 



dV L 
> ~df 



dV L 
dt 



i{[H s + H Ll V L ]) =0. (7) 



Note that we have assumed that [V L , Vr] — 0. Using Eq. 
([7]) we reach the following expression for the averaged 
energy current (valid in the long time limit), 

(J L ) = Tr s Tr F [J L p(i)] - -i([H s , V L \). (8) 

This commutator can be readily evaluated to yield, 

[Hs,Vl] = 2iAa y ^^2aL,i:LMC Ll CL,i', (9) 



1,1' 



leading to 



(J L ) = 2ATr s [a y Tr F [A L p(t)}} 



(10) 
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with the bath operator Al = J2i v a L,i-L,i'c\ iCl,v- We 
now define a subsystem operator 

A s (t) = Tr F [A L p(t)] = Tr F [e iHt A L e- im p(0)], (11) 

and express the current using its matrix elements 

(J L ) = 2A[-i(4sr(t))i,o + i(4s(i))o,i]. (12) 

We emphasize that this expression is designed to provide 
the steady-state value and not the transients, given our 
assumption ([7]). 

The two operators, ps(t) and As{t), are subsystem op- 
erators. They are simulated in the next section directly, 
using INFPI, or studied in a perturbative manner, under 
the Markovian limit, to provide Kinetic-type expressions. 

IV. METHODS 
A. Path-integral simulations 

The principles of the INFPI approach have been de- 
tailed in Refs. [l(l[l7]]> where it has been adopted for in- 
vestigating dissipation effects in the nonequilibrium spin- 
fermion model and charge occupation dynamics in the 
interacting Anderson model. Other applications include 
the study of the intrinsic coherence dynamics in a double 
quantum dot Aharonov-Bohm interferometer [3lj , explo- 
ration of relaxation and equilibration dynamics in finite 
metal grains [32j |. and the study of electron-phonon ef- 
fects in molecular rectifiers [331 ] . 

Here, using INFPI, we can directly simulate both the 
dynamics of the reduced density matrix ps{t) [Eq. J5J], 
and the time evolution of subsystem expectation values, 
particularly As(t) [Eq. (|TTj) ]. which can be used to ob- 
tain the energy current (Jl), Eq. p^|) . In practice, for 
achieving fast convergence, we have simulated directly 
the averaged current 

J = \{{Jl) - {Jr})- (13) 

The negative sign in front of (Jr) arises from our sign 
convention; the current (J„) is defined from the v reser- 
voir, into the junction. 

Algorithmic details of the INFPI method were recently 
elaborated in Ref. [33], thus we only include the main 
principles here. The algorithm is based on a Trotter 
breakup of a short-time time evolution operator into two 
parts: a (simple) time evolution term that depends on 
the subsystem Hamiltonian, and a term that accommo- 
dates the reservoirs Hamiltonians, and their interactions 
with the subsystem. Collecting the contribution of the 
latter terms along the subsystem path, we construct the 
so called "influence functional" (IF), which involves non- 
local dynamical correlations. The IF has an analytical 
form in some special cases jisj ; in the present model its 
form is only known in the weak-intermediate coupling 
limit [H, Q , thus we evaluate it numerically by energy- 
discretizing the Fermi sea. 



The main conceptual element behind the INFPI ap- 
proach is the observation that at finite temperatures 
and/or nonzero chemical potential bias bath correlations 
exponentially decay in time, allowing for their trunca- 
tion beyond a memory time t c . The dynamics can then 
be achieved by defining an auxiliary density matrix, or 
more generally, a subsystem operator [e.g., As(t) of Eq. 
dill) ], on the time- window t c . This nonlocal object can 
be iteratively evolved from the subsystem-bath factorized 
initial condition, to the present time t. 

This path-integral method involves three numerical pa- 
rameters: (i) the number of states used in the discretiza- 
tion of each Fermi sea L, (ii) the time step adopted in the 
Trotter breakup St, and (iii) the memory time accounted 
for, t c , beyond which the IF, accommodating the effect 
of the reservoirs on the subsystem, is truncated. Conver- 
gence of INFPI is verified by confirming that results are 
insensitive to the reservoirs discretization, the finiteness 
of the time step, and the memory size t c = N s St, with 
N s as an integer. It should be noted that minimizing 
the Trotter breakup error, taking St — > 0, conflicts with 
the need to cover the memory time- window t c , since the 
parameter N s has to be increased until convergence is 
reached. Since our computational effort scales as d 2Ns , 
where d is the Hilbert space dimensionality of the sub- 
system, we are practically limited to N s < 10, implying 
on the minimal time step that can be adopted. 



B. Markovian Master Equation 

The dynamics of the model @, and its variants, can 
be analyzed in the weak subsystem-bath coupling limit 
under the Markovian approximation [l3|, HH . The prob- 
abilities P n to occupy the |n) state of the subsystem, the 
quantum impurity, satisfy the Master equation 

Pn ^ Prnkjn^-n Pn ^ ^n— >m> (1^0 

m m 

where the transition rate from the state |to) to \n) (to ^ n 
and to, n—0,1 here) is additive in the L and R reservoirs, 
k n ^ m = + kn-+ m , due to the linear form of the 

interaction |2lL l3o| . 

This type of Kinetic equation has been used in many 
recent works for investigating energy, spin, and charge 
transfer in open quantum systems [2Ci l2l| . Particularly, 
it has been recently adopted for modeling radiative en- 
ergy transfer between metals [13l |28| , and for studying 
charge and energy transfer phenomenology in mesoscopic 
systems [lj| HH [H| and single molecules [23[ . It is 
thus important to test the suitability and accuracy of this 
common and well-accepted approximate scheme against 
exact results. 

In steady-state, taking Eq. (|12p as a starting point, 
one can show that in the weak coupling limit and under 
the Markovian approximation the energy current across 



4 



the system reduces to 29] ((Jl) = J in steady-state), 

J = ^ E mtTl P n k^ m , (15) 



with E„ 



E m — E n . The current is defined positive 



when flowing left to right. At the level of the Golden- 
Rule formula, the transition rates are given by [l~3T ] 

h v = 

2Tr^2\a v ,j; Vtj >\ 2 n u F (e k )[l - n" F {ej>)]5(tj - e r - E m<n ) 

3 J' 



2tt J den F (e)[l-n F (e-E m , n )]F u (e). 



= -2im B (E m ,„) J de [n F (e) -n F (e- E m , n )] F„(e)(16) 

From the last relation we note that the thermal properties 
of the reservoirs are concealed within both the Fermi- 
Dirac distribution function n F {e) — [e( e ~ M, ^/ Tl ' + 
and the Bose-Einstein occupation factor n B (e) = [e e / T " — 
l] -1 . It is therefore clear that when the integral yields 
a temperature independent constant, the statistic of the 
reservoirs is fully bosonic [Hj]. The other element in Eq. 
([TBI) is a dimensionless interaction term 



F v (e) = \a v \ 2 p v {e)p v (e - E m ,n), 



(17) 



which encloses the properties of the reservoirs, multiplied 
by the subsystem-bath (energy independent) couplings 
a v . Once we assume that the density of states is a con- 
stant 36], F u (e) « F v (fi v ), the integration in Eq. (fTF|) 
can be performed when the Fermi e nerg ies are situated 
far from the conduction band edges [36J. Making use of 
the following relation, 



de[n F (e) - n F (e - F m ,„)] = -E m>n , (18) 



we reach the closed-form expression, 

K^ m = 2nn B {E m . n )E m ^ n F u (p, u ). (19) 

Note that n B {—E m ^ n ) = ~[n B (E m ^ n ) + 1], thus the ex- 
citation and relaxation rates induced by the v reser- 
voir satisfy the detailed balance relation, ^^ rm jk^ n ^ tn = 
e -s m _„/T„_ can a j g0 ex p ress th e ra tes in terms of a 

subsystem-bath interaction parameter 



T£(2A) = 2ttF,(^)2A 



2A o 
2 [Trp^K] 2 , (20) 

IT 



where we recall that both the density of states and the 
interaction parameter a are assumed to be energy inde- 
pendent. Using this defining, the rate constants in our 
model reduce to 



K l- 

" 11. 



0^1 



r£(2A)[l + ra£(2A)], 
r^(2A)n£(2A). 



(21) 



For simplicity, we do not include below the explicit de- 
pendence of T F and n B on energy; both quantities should 



be evaluated at the subsystem energy gap 2 A. We cal- 
culate the population of the states in steady-state by 
putting P n = in Eq. (Q3]). With this at hand, the 
energy current (fT5)) simplifies to [37J 



J = 2A 



nn K 



r£(l+2n£)+r«(l + 2ng) 



(22) 



This expression provides the steady-state energy current 
in the weak coupling limit, under the Markovian approx- 
imation. 

It should be noted that a more involved non-interacting 
blip-approximation (NIBA) type scheme [l[ can be imple- 
mented for following the qubit dynamics in the nonequi- 
librium spin-fermion model [7J. Furthermore, besides the 
qubit dynamics itself, the heat current can be simulated 
within the NIBA approximation by extending a gener- 
ating function technique developed in Ref. [38j for the 
study of transport behavior in the nonequilibrium spin- 
boson model. However, here we contain ourselves with 
the simpler, more standard and common Markovian Mas- 
ter equation, with the objective to provide insight on 
its applicability and accuracy for the large community 
adopting it in studies of charge, spin, and energy dynam- 
ics in mesoscopic and molecular systems. 
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FIG. 2. Energy current as a function of metal-qubit coupling 
parameter Tp using the bandwidth D — 2, A = 0.1, Tl = 0.4, 
and Tr — 0.2. INFPI numerical parameters are St = 1 and 
7V S = 9. Dashed line: Master equation results. INFPI results 
appear in symbols, □ for L — 40 and o when taking the 
asymptotic L — > oo limit. Inset: Zooming over the small 
coupling limit where the current linearly scales with T F . 



V. RESULTS 

We compare INFPI simulations to the Master equa- 
tion predictions. Within INFPI, the current is simulated 
directly using Eqs. (fTTj) and (fT2"j). and we show results 
only in the long time (quasi) steady-state limit. The 
closed- form Master equation expression is given by (|22[) . 
Beyond the weak coupling limit, we calculate the Master 
equation current directly from Eq. (| 15[) using the rates 
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FIG. 3. Converging the data of Fig. Uto the ->■ limit, 
with the intercept representing the asymptotic result. 
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FIG. 4. Polarization dynamics at different coupling strengths, 
for the same set of parameter as in Fig. [2] Results are dis- 
played in the basis of the Hamiltonian Q 



(|16[) . since the wide-band approximation does not hold 
anymore (though of course the applicability of the Mas- 
ter equation technique as a whole is questionable in this 
regime) . 

We typically use the following set of parameters: an 
energy gap 2A = 0.2 (arbitrary units) for the subsystem, 
metallic bandwidth D = 2, T v > A, and the equilibrium 
Fermi energies located at the center of the band. We also 
define the dimensionless parameter 



% = irp u a Ul 



(23) 



which is varied between and 0.8 here, where conver- 
gence is achieved. This choice corresponds to < r^, < 
0.08 when A = 0.1, see Eq. pp]). For this set of model pa- 
rameters, we have confirmed that selecting St = 0.8 — 1.5 
and 7V S = 6 — 9 (yielding memory time r c > 8) provides 
converging results, see panel b in Figs. El Eland [3 

Before turning to simulations, it is important to iden- 
tify the region of weak subsystem-bath coupling. It holds 



when the scattering phase shift is small, arctan<^„ ~ <j) v 
39]. This translates (within 5% error) to the dimension- 
less parameter <j> being limited to the domain of </>„ < 0.4. 
Equivalently, the interaction energy should be limited to 
T F < 0.02, or in other words, T F /(2A) < 0.1, i.e., in the 
weak coupling limit the subsystem gap is large relative 
to the coupling energy. 

Fig. [2] displays the energy current as a function of the 
interaction energy for a symmetric junction, T F = T F . 
We find that when T F > 0.02, Master equation-derived 
energy current overestimates the exact result by more 
than 10%. In the strong coupling limit, Master equa- 
tion overvalues the correct numbers by a factor of two. 
More significantly, this Golden-Rule based method can- 
not reproduce the saturation effect of the current with the 
subsystem-bath interaction parameter, and it wrongly 
predicts a linear scaling, J oc T F . Note that we include 
INFPI results both for the case of L = 40 bath states (□), 
and in the asymptotic L — ¥ oo limit (o), obtained by ex- 
trapolating the linear J vs. L _1 curves to L^ 1 0, see 
Fig. [3l We find that this extrapolation affects the results 
by up to 4% at strong coupling, while the weak coupling 
values are unaffected. While we do not show transient 
data for the current, we comment that steady-state has 
been reached at At ~ 50 — 100 in the weak coupling 
limit; it is established much faster, At ~ 5 — 10 at strong 
coupling. 

We correlate transport behavior of the junction with a 
study of the dissipative dynamics of the spin polarization, 
as obtained from INFPI, in Fig. SJ We observe weakly- 
damped coherent oscillations in the weak coupling limit 
when perturbative Master equation well describes the dy- 
namics. These oscillations still survive at intermediate 
coupling, but at strong coupling the polarization expo- 
nentially decays in time (inset). Crucial parameters of 
the model are the scattering phase shifts S±. In equilib- 
rium, the phase shifts are given by [1^, |4Cj 



tan<5± = 4>l,r 



(24) 



Out-of-equilibrium, A/x ^ 0, the phase shifts are com- 
plex numbers [4oj ]. In the spin-boson model the Kondo 
dimensionless dissipation coefficient £ represents a char- 
acteristic exponent in the system: at zero temperature 
and zero energy bias the spin displays damped coher- 
ent oscillations for £ < 0.5, relaxation dynamics between 
0.5 < £ < 1, and a localization phase for £ > 1 [4l| . Thus, 
this parameter controls dissipation-induced phase transi- 
tions. In the fermionic analogue it can be shown that 
the characteristic exponent is given by £= (5'1 + 5 2 _)/tt 2 
[H- Since \S±\ < tt/2, £ < 1/2. Thus, in the spin- 
fcrmion model described in this paper the spin cannot 
manifest the localization behavior, and a large <f> value 
brings us to the relaxation scenario, as indeed observed 
in Fig. U] 

The current-temperature characteristics of the junc- 
tion are depicted in Figs. EH using different coupling 
strengths. We find that at weak coupling Markovian 
Master equation very well reproduces the qualitative and 
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FIG. 5. (a) Energy current-temperature characteristics for 
the same set of parameter as in Fig. [2] (j> = 0.2 (r^ = 0.005). 
We vary Tl, but keep Tr fixed, Tr = 0.2. (b) Convergence 
behavior with increasing memory size. Data was produced 
with three different time steps, St = 0.8, f.O and 1.2. 
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FIG. 6. (a) Energy current-temperature characteristics for 
the same set of parameter as in Fig. [2] cj> = 0.5 (Tp = 0.032). 
We vary Tl, but keep Tr fixed, Tr = 0.2. (b) Convergence 
behavior with increasing memory size. Data was produced 
with three different time steps, St = 0.8, 1.0 and 1.2. 



quantitative aspects of the current, even at a large tem- 
perature difference, see Fig. [5] The convergence behav- 
ior of the INFPI method at different temperature biases 
is displayed in panel (b) , where we show the energy cur- 
rent as obtained using different memory time r c and time 
steps. At intermediate couplings, Fig. |6]shows that Mas- 
ter equation overestimates exact results by up to 25% at 
large temperature differences. When the subsystem-bath 
coupling is large, we managed to converge simulations 
only up to the bias T L — Tr ~ 0.2. The Kinetic method 
now provides values that are a factor of 2 larger than the 
exact numerical data. It is important to note that the 
qualitative current-temperature features are correctly re- 



FIG. 7. (a) Energy current-temperature characteristics for 
the same set of parameter as in Fig. [2] <j> — 0.75 (Ff = 0.072). 
We vary Tl, but keep Tr fixed, Tr = 0.2. (b) Convergence 
behavior with increasing memory size: results converge only 
at small bias, Tl — Tr < 0.4. Data was produced with three 
different time steps, St = 0.8, 1.0 and 1.2. 



produced within the Markovian Master equation, even at 
strong coupling. However, if one is interested in quan- 
titative information, Master equation can be used in its 
strict regime of applicability only, <fi < 0.2. Data in Figs. 
[5][7] is presented for a fixed value of bath states, L = 40, 
since we have confirmed that taking the large-L limit 
only corrects the current by < 4%, at both low and high 
temperature biases. 



VI. SUMMARY 

We have studied energy transfer between metals me- 
diated by a quantum impurity, using two approaches: 
numerically exact path-integral simulations and analytic 
results from a Golden- Rule type Markovian Master equa- 
tion treatment. We found that standard Master equa- 
tions fail to reproduce the current-interaction energy 
characteristics already at intermediate system-bath cou- 
plings, as it can only provide a linear enhancement of 
the current with the subsystem-bath interaction, missing 
a saturation effect. In contrast, the current-temperature 
characteristics is produced in a qualitative correct way 
by a Master equation formalism, though actual values 
deviate by 100%, and more, at high temperature biases 
and at strong coupling. 

Our results are beneficial for the critical testing of com- 
mon Master equation techniques. The methods described 
are also useful for practically modeling superconducting- 
based qubit devices [42| ■ While a Master equation treat- 
ment offers simple-intuitive expressions that often allow 
to discern essential transport characteristics, already at 
intermediate system-bath couplings it may overestimate 
the current by ~ 10%, up to a 100% incorrect enhance- 



7 



ment at strong coupling. These deviations are certainly 
important when a quantitative analysis of device effi- 
ciency is performed. In particular, the calculation of en- 
ergy conversion efficiency in conducting junctions should 
be done with caution when a Master equation is of use 

[mum. 

In our future work we plan to study the heat current 
characteristics in the complementary spin-boson type 
molecular junction model [33]. This could be done by 
extending the Feynman- Vernon IF expression [44[ to de- 
scribe the evolution of other operators besides the re- 



duced density matrix. Alternatively, one could use the 
bosonization approach and draw general results for the 
nonequilibrium spin-boson problem, based on the spin- 
fcrmion model calculations presented here. 
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